#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdbool.h>
#include <quadmath.h>

#include "../src/myexp.h"
#include "../src/myexp_impl.h"

#define ARY_SIZE(a) sizeof(a)/sizeof(a[0])

typedef void (*fp_v4_exp)(double *pr, const double *px );

void print_result_source_ary_len4(const double xs[4],const double ys_ref[4],const double ys[4])
{
  bool r=true;
  for (int i=0;i<4;i++)
  {
      printf("exp(%lf)=%.13la:exp_ref\n",xs[i],ys_ref[i]);
      printf("exp(%lf)=%.13la:exp_vec4\n",xs[i],ys[i]);
  }
  printf("\n");
}

void test_expvec4_case0(fp_v4_exp fp)
{
  printf("\nrun %s()\n", __func__);
  double xs[4]={
    (1+0.0/128.0)*LN2,
    (1+1.0/128.0)*LN2,
    (1+2.0/128.0)*LN2,
    (1+3.0/128.0)*LN2
  };

  double ys[4];
  double ys_ref[4];

  for (int i=0;i<4;i++)
  {
    double x = xs[i];
    ys_ref[i] = exp(x);
  }
  fp( ys, xs );
  print_result_source_ary_len4(xs,ys_ref,ys);
}

void test_expvec4_case1( fp_v4_exp fp)
{
  printf("\nrun %s()\n", __func__);
  double xs[] =
  {
    -1.0/0.0 , 1.0/0.0 /* nan*/,
    -__builtin_inff(), __builtin_inff(),
  };

  double xs2[] =
  {
    -1025.0, -1024.0, -1023.0, -1022.0, -1021.0,
     1021.0,  1022.0, 1023.0,  1024.0,   1025.0
  };

  double block[4];
  double ys_ref[4];
  double ys[4];

  double *ary;
  double scale;
  int ary_size;
  for (int round=0;round<2;round++)
  {
    ary_size=(round==0)? ARY_SIZE(xs): ARY_SIZE(xs2);
    ary=(round==0)? xs: xs2;
    scale=(round==0)? 1.0: LN2;

    for (int i=0;i<ary_size;i++)
    {
      for (int j=0;j<4;j++)
      {
        block[j] = (j==(i & 3)) ? ary[i]*scale : 1.0;
        ys_ref[j]= exp(block[j]);
      }
      fp( ys, block );
      print_result_source_ary_len4(block,ys_ref,ys);
    }
  }

  ary_size=ARY_SIZE(xs2);
  ary=xs2;
  scale=LN2;
  for (int i=0;i<ary_size;i++)
  {
    for (int j=0;j<4;j++)
    {
      block[j] = (j==(i & 3)) ? (ary[i]+0.001)*scale : 1.0;
      ys_ref[j]= exp(block[j]);
    }
    fp( ys, block );
    print_result_source_ary_len4(block,ys_ref,ys);
  }

}

void test_expvec4_case2(fp_v4_exp fp)
{
  printf("\nrun %s()\n", __func__);
  int e;
  double f,eps,x,x1,y0,y1;

  double xs[4];
  double ys[4];
  double ys_ref[4];

  //1022*log(2) ~=  708.39641853226410
  for (int i=0;i<1000;i++)
  {
    for (int j=0;j<4;j++)
    {
      int x= rand() % 141600;
      x -= (141600/2);
      xs[j]= (double)(x * 0.01); // fabs(x*0.01)<708.00
      ys_ref[j] = exp( xs[j] );
    }
    fp(ys,xs );
    print_result_source_ary_len4(xs,ys_ref,ys);
  }
}

void test_expvec4_case3(fp_v4_exp fp)
{
  printf("\nrun %s()\n", __func__);
  double f,x,x1,y0,y1;

  double xs[4];
  double x1s[4];
  double ys[4];
  double ys_ref[4];

  for (f=-10.0; f<=10; f= f + (0.01*4) )
  {
    for (int j=0;j<4;j++)
    {
      x  = f + (j * 0.01);
      x1 = x * INV_LN2;
      xs[j]=x;
      x1s[j]=x1;
      ys_ref[j] =(double)expq((__float128)x);
    }
    fp(ys,xs );
    print_result_source_ary_len4(xs,ys_ref,ys);
  }
}

void test_exp_vec4_avx_fun()
{
  printf("Run in %s\n",__func__);
  test_expvec4_case0(exp_vec4_avx);
  test_expvec4_case1(exp_vec4_avx);
  test_expvec4_case2(exp_vec4_avx);
  test_expvec4_case3(exp_vec4_avx);
}

void test_exp_vec4_avx2_fun()
{
  printf("Run in %s\n",__func__);
  test_expvec4_case0(exp_vec4_avx2);
  test_expvec4_case1(exp_vec4_avx2);
  test_expvec4_case2(exp_vec4_avx2);
  test_expvec4_case3(exp_vec4_avx2);
}

int main(int argc, char *argv[])
{

#ifdef _DEBUG_
  //debug_exp();
  return 0;
#endif
  test_exp_vec4_avx_fun();
  
  if ( __builtin_cpu_supports("fma") ) 
  {
      test_exp_vec4_avx2_fun();
  }
  return 0;
}
